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Abstract 

We present an algorithm that computes the geodesic center of a given polygonal domain. 
The running time of our algorithm is for any e > 0, where n is the number of corners 

of the input polygonal domain. Prior to our work, only the very special case where a simple 
polygon is given as input has been intensively studied in the 1980s, and an 0(n log n)-time 
algorithm is known by Pollack et al. Our algorithm is the first one that can handle general 
polygonal domains having one or more polygonal holes. 


1 Introduction 

The diameter and radius of a compact shape are among the most natural and fundamental param¬ 
eters describing and summarizing the shape itself. In this paper, we study these quantities for a 
polygonal domain V, that is, a polygon having h>t) holes. More specifically, a polygonal domain 
is a connected and compact subset of whose boundary consists of /i + 1 simple closed polygonal 
curves. In regard to a metric d on. V, the diameter of V is defined to be the maximum distance over 
all pairs of points in the V, that is, maxp^^g-p d{p, q), while the radius is defined to be the min-max 
value minpgp max^gp d{p, q). A pair of points in V realizing the diameter is called a diametral pair, 
and a center is defined to be a point c G U such that maxggp d{c, q) is equal to the radius. Among 
common metrics on a polygonal domain V, we consider the geodesic distance d(p, q) for p,q £ V 
that measures the Euclidean length of a shortest path that connects p and q and stays inside V. 
The diameter and the radius of a polygonal domain V with respect to the geodesic distance d are 
often called geodesic diameter and geodesie radius of V, respectively. 

The problem of computing the geodesic diameter and radius of a simple polygon (i.e., a polyg¬ 
onal domain with no holes) has been intensively studied in computational geometry since the early 
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80s. For the geodesic diameter problem, Chazelle gave the first algorithm whose running time 
was O(n^), where n denotes the number of vertices or corner^ of the input polygon. This was 
afterwards improved to 0(n log n) time by Suri [^, and hnally to linear time by Hershberger and 
Suri [^. For the geodesic radius of a simple polygon, the first algorithm was given by Asano and 
Toussaint [^, and its running time was 0(n^ log n)-time. Later Pollack, Sharir, and Rote 
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proved it to 0(nlogn) time. Very recently, an optimal 0(n)-time algorithm for the geodesic radius 
of a simple polygon is presented by Ahn et al. [^. 

The case in which the domain has one or more holes is much less understood. To the best of our 
knowledge, the only known result is a companion paper in which an algorithm that computes the 
geodesic diameter of a polygonal domain with n corners and h holes in or 0(n'^(log n + h)) 

time is given. As for computing the radius, no algorithm was known prior to this work, even 
though the problem has been remarked repeatedly as an important open problem 11, Open Problem 
6 ]. 


The main difference between simple polygons and general domains lies on the difficulty to 
determine and discretize the search space. A key tool often used in these problems is, given a point 
p € V, compute a farthest neighbor of p, a point of V that is farthest away from p. It is well known 
that in simple polygons, every farthest neighbor of any point V should be a corner of P [^. This 
implies that the geodesic diameter of any simple polygon can only be determined by two of its 
corners. In particular, the problem is now reduced on how to efficiently search among the O(n^) 
candidates that can potentially determine the geodesic diameter, so one could try any bruteforce 
search on them. The geodesic radius of a simple polygon can also be handled in a similar way: 
Even though the corresponding center itself may be an interior point of V, its farthest neighbors 
are all corners. 

For general polygonal domains, unfortunately, this is not the case any more. A farthest neighbor 
of a point in a polygonal domain V having one or more holes may not be a corner of V, and even can 
be an interior point of V. This makes things complicated; the geodesic diameter can be determined 
by two interior points, as shown in [^. This difference mainly causes the huge gap, 0(n) and 
in the computational complexity of computing the geodesic diameter between simple 
polygons and general domains [^. 

In this paper, we present an algorithm that, in time, computes the geodesic radius 

and center. At a glance, the time complexity might seem very high, but it is comparable to 
the currently best algorithms for computing the geodesic diameter. Indeed, a crucial observation 
for the diameter algorithm is that there are at least hve shortest paths between the two points 
determining the geodesic diameter if the two points lie in the interior of R [^. This observation 
leads to a bounded number of candidates for diametral pairs. In Section we show that the 
geodesic radius and center sometimes involves nine paths to determine. This enlarges the search 
space considerably, thus a larger running time is somehow expected. 

The rest of the paper is organized as follows. After introducing preliminary definitions and 
concepts in Section we list geometric observations in Section that will be the base of our 
algorithm described in Section Finally, Section concludes the paper with possible lines of 
future research. 


corner of a polygon usually indicates a vertex to which two incident edges form an angle that is not 180° 
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(a) (b) (c) 

Figure 1: (a) A polygonal domain V with two holes and a source point s £ V. (b) The shortest 
path tree SPT(s) on V U{s} with root s (edges are directed towards descendants), (c) The shortest 
path map SPM(s) (depicted by solid segments). Corners v £ V with non-empty region <Ts(n) are 
marked by black dots. 


2 Preliminaries 


Throughout the paper, we frequently use several topological concepts such as open and closed 
subsets, neighborhoods, and the boundary dA of a set A] unless stated otherwise, all of them are 
derived with respect to the standard topology on with the Euclidean norm || • || for fixed d> 1. 
We also denote the straight line segment joining two points a,b £M? by ab. 

A polygonal domain V with h holes and n corner^ is a connected and closed subset of with 
h pairwise disjoint holes. Each hole is a simple polygon contained in V. Thus, the boundary dV of 
V consists of h + 1 simple closed polygonal chains, and overall n line segments. Each of the holes 
(and the outer boundary ofV) is regarded as an obstacle that feasible paths in V are not allowed to 
cross. The geodesic distance d{p, q) between any two points p,qfn. a polygonal domain V is defined 
to be the Euclidean length of a shortest feasible path between p and g, where the length of a path 
is the sum of the Euclidean lengths of its segments. It is well known 10 that there always exists 


a shortest feasible path between any two points p,q £V, and the geodesic distance function d(-, •) 
is thus well defined. 

The geodesic radius rad('P) of V is defined to be the min-max quantity: 


rad('P) = minmaxd(p, g). 

p&V q&V 


A geodesic center of "P is a point c£V such that 


maxd(c, q) = rad(P). 
q£P 

The set of all geodesic centers of V is denoted by cen(P). The purpose of this paper is to 
describe the first algorithm that exactly computes the geodesic radius rad(P) and centers cen(P) 
of a given polygonal domain V. 


2.1 Shortest path trees and shortest path maps 

Let V be the set of all corners of P and tt be a shortest path between any two points s,t £ P. 
This path vr is a polygonal chain that makes turns only at corners V of P [^. We represent tt 

^We reserve the term “vertex” for a 0-dimensional face of subdivisions of a certain space. 
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by the sequence of traversed corners: tt = (s, ui,..., t) for some vi,... ,Vk G V. Note that k 
may be zero; in this case, the shortest path tt is the segment st connecting the two endpoints, and 
thus d{s,t) = ||s — 1 ||. If two paths (with possibly different endpoints) induce the same sequence of 
corners (ui ,... ,Vk), then they are said to have the same combinatorial structure. 

Given a source point s G "P, the shortest path tree SPT(s) of s is a tree spanning V U {s} 
embedded in V such that the unique path in SPT(s) from s to each corner of P is a shortest path 
in P. See for example Figure [^b). 

The shortest path map SPM(s) of a fixed s G P is a decomposition of P into cells such that 
points in the same cell can be reached from s by shortest paths of the same combinatorial structure. 
See Figure[2c). Each cell as{v) of SPM(s) is associated with a corner u G E which is the last corner 
of the shortest path tt from s to any t in the cell crs(v). Note that the path vr goes along the path 
in SPT(s) to V and then reaches t along vt. We also define the cell crs(s) as the set of points t G P 
such that st C P, i.e., d(s, t) = ||s — t\\. 

Edges of SPM(s) either belong to &P or are arcs on the boundary of two incident cells crs(ui) 
and as{v 2 ) determined by two corners ui,U 2 G E U {s}. Edges of the second kind are hyperbolic 
arcs if ui and V2 are not adjacent in SPT(s). Moreover, there are two different shortest paths from 
s to any point on an edge of SPM(s), one via vi and the other via V2. 

Vertices of SPM(s) are either corners of P, endpoints of an edge of the second kind above, or a 
point p G V incident to at least three faces cts( ui), (Ts(u 2 ), (Ts(u 3 ) for some corners vi,V 2 , U 3 G EU{s}, 
yielding three different shortest paths from s. Depending on which of the three cases it falls into, 
each vertex of SPM(s) admits either 1 , 2 , or more different shortest paths from s to the vertex, 
respectively. 

The shortest path map SPM(s) has 0(n) cells, edges, and vertices in total, and can be computed 
in 0(nlogn) time using 0(n logn) working space [^. For more details on shortest path maps, we 
refer to [9[|Tl] . 

2.2 Path-length functions 

For any point p G P, we define its visibility region as the set VR(p) of all points q G V such that 
pq CV, that is, points q that sees p. 

Let vr be a shortest path from s to t for s,t gV. If tt 7 ^ st, then there are two corners u,v G V 
such that u and v are the first and last corners along vr from s to t, respectively. Here, the path 
TT is formed to be the union of Mt, vt and a shortest path from u to v. Note that u and v are 
not necessarily distinct. In order to realize such a path, s must see u and t must see v, that is, 
s G VR(u) and t G VR(u), 

We define the path-length function len„^^: VR(u) x VR(u) —>■ M for any fixed pair of corners 
tt, u G E to be 

len„^t,(s, t) := ||s — u|| + d(tt, v) + ||u — t\\. 

That is, lenu^^(s,f) represents the length of paths from s to f that have a common combinatorial 
structure; going straight from s to tt, following a shortest path from tt to u, and going straight 
to t. Also, unless d{s,t) = ||s — t\\ (equivalently, s G VR(t)), the geodesic distance d{s,t) can be 
expressed as the pointwise minimum of path-length functions. 

d(s,t) = min lenu„(s,t). 

«eVR(s)ny, v&VR{t)nv 
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By definition of shortest path map SPM(s) and its cells crs{v), if t G for some v £ V, then 
we have d{s,t) = lenu^^(s,t), where u £V denotes the first corner along the shortest path from s 
to V, or equivalently, along the path from s to u in SPT(s). 


3 Farthest Neighbors and Geodesic Centers 


In this section we introduce several tools that will be useful for discretizing the search space. For 
any point p £ V, we let <h(p) be the maximum geodesic distance we can obtain when we fix one 
point as p, that is, 

<I>(p) := maxd(p, g). 
qeV 

We call a point q £ V a farthest neighbor of p G "P if d{p, q) = <h(p). 

Observe that the geodesic radius of V is the minimum possible value of <h(p) over all p £ V, 
that is, 

rad(P) = min^>(p), 
pS'P 

and a point that minimizes $(;>) is a geodesic center of V. 

The following lemma gives us a way of computing farthest neighbors. Recall that each vertex 
of the shortest path map SPM(p) for p £ V is either a corner of P, an endpoint of an edge lying on 
the boundary dV, or a point in the interior of P that is incident to three edges. 


Lemma 1. For any point p £ V, any farthest neighbor of p in V is a vertex o/ SPM(p). 


Proof. Suppose for the sake of contradiction that there exists a farthest neighbor q £V oi p that 
is not a vertex of SPM(p). Then, there exists a sufficiently short line segment L such that L is 
contained in the closure of some cell (Tp{v) of SPM(p) for some u G R U {p} and contains q in its 
relative interior. This is always true even if q lies on an edge of SPM(p) since every edge of the 
shortest path map is either straight or hyperbolic. 

Then, the function f{x) = d{p,x) for x G L is represented as f{x) = len„^„(p, x) = \\p — u|| + 
d(u,u) + ||x — x|| for some u £ V L) {p}. Observe that the function / is convex on L and has no 
plateau along its graph. Since q lies in the relative interior of L, there always exists a point y £ L 
such that f{y) > f{q), which contradicts the assumption that g is a farthest neighbor of p. 0 


This observation is analogous to the fact that farthest neighbors of any point in a simple polygon 
However, vertices of shortest path maps SPM(p) may lie in the interior of P. 
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are its corners 

This means that a geodesic radius and center may be determined by interior points, whereas this 
never happens for simple polygons. 

Figure illustrates an example polygonal domain such that its unique geodesic center c and 
its farthest neighbors lie in its interior. The domain shown in Figure consists of three identical 
regions arranged in a symmetric way: each part contains two holes that almost ht together forming 
a very narrow corridor between them. We claim that c is the unique geodesic center and c has 
exactly three furthest neighbors: qi, q 2 , and q^. 

Observe hrst that each q^ is a vertex of the shortest path map SPM(c) for c. Each vertex of 
SPM(c) lying in the interior admits three distinct shortest paths from c. (Moreover, because of 
the way the regions are defined, each q^ is also the farthest neighbor of c within the corresponding 
region. Due to symmetry in the construction, we observe that no point other than c can be closer 
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Figure 2: A polygonal domain instance with a unique interior geodesic center c. The three farthest 
neighbors qi, q 2 , and q^ of the center c lie in the interior of the domain as well. Observe that there 
are three distinct shortest paths between the center c and each of its farthest neighbors, and thus 
9 shortest paths of equal length in total. 


to all of qi, q 2 , and gs at the same time. Thus, the point c is the only geodesic center of this 
polygonal domain. 

Note that this construction is slightly degenerate since it has a few symmetries. However, such 
degeneracies can be removed by a small perturbation on the location of the corners. This is possible 
because the center c is “stable” in the sense that a sufficiently small perturbation on the corners 
of the domain will only imply a small change in the location of c and points g*. 

4 Algorithm 

In this section, we describe our algorithm for computing the radius rad('P) and all centers cen('P) 
of an input polygonal domain V. Recall that the problem of computing the radius and center can 
be seen as a minimization problem under the objective function $ over V. Thus, our approach is 
to decompose V into cells, and find candidate centers in each cell. 

For any subset u C P of the domain V, we call the minimum value of ‘h(p) over p £ a 
the a-constrained geodesic radius, and each point in a that attains the minimum is a a-constrained 
geodesic center. Clearly, in any decomposition {cJi, cr 2 ,...} of "P, the geodesic radius is the minimum 
of CTj-constrained geodesic radii over all i, and the points that attain the minimum value form the 
geodesic centers cen(P) of P. 

In this paper, we use the SPM-equivalence decomposition. This decomposition subdivides a 
polygonal domain P into cells such that for all points s in a common cell a of ^spMj their shortest 
path maps SPM(s) are topologically equivalent. More precisely, two shortest path maps SPM(si) 
and SPM(s 2 ) are said to be topologically equivalent if their underlying labeled plane graphs are 
isomorphic. This structure was introduced by Chiang and Mitchell as a means to devise efficient 
data structures that support two-point queries for Euclidean shortest paths. In their work, they 
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show that the decomposition ^spm has complexity and can be computed in 0(n^° logn) 

time. 

An additional property of this subdivision (also shown by Chiang and Mitchell (^) is that, for 
any cell a of AspM) the elements of SPM(s) (i.e., vertices and edges) can be explicitly described by 
algebraic functions of s G u. In this manner, the shortest path map SPM(s) for any point s G cr 
can be parameterized within a hxed cell a. 

Let a be any cell of Aspm- Since all shortest path maps SPM(s) within s G cr are topologically 
equivalent, they must all have the same number m of vertices. We are particularly interested in 
coordinates of the vertices vi,V 2 , ■ ■ ■ ,Vm of SPM(s) as functions of s G cr. Recall that the vertices 
Vi of SPM(s) must include the corners V of V] thus, if Uj G R for z = 1, 2,..., m, then Vi{s) will be 
a constant function that maps to a unique corner of V. 

For i = 1, 2,..., m, we define the function /j: cr —>• M to be fi{s) = d(s, Ui(s)) for s G a. That is, 
this function maps s to the geodesic distance from s to Vi{s). We then consider the upper envelope 
maxj fi{s) of the m functions, which maps s to its maximum geodesic distance over all the vertices 
Vi{s) of SPM(s). By Lemmathe farthest neighbors of s must be among the Vi{s), and it thus 
holds that 

$(s) = max fi{s). 

2=1,2,...,m 

In order to find the cr-constrained geodesic radius and center, it suffices to compute and search the 
upper envelope of the m functions /*. 

In order to obtain an explicit expression of fi{s), we observe the following property. 

Lemma 2. For any cell a of Aspm o-nd i G {1, 2,..., m}, one of the following holds for all s G a: 
a C VR(ui(s)) or (T n VR(ui(s)) = 0. 

Proof. This follows from the fact that shortest paths are topologically equivalent within a. If 
s G a sees Vi{s), then the shortest path tt from s to Vi{s) is sui(s). Since the shortest path vr' from 
any other s' G cr to Uj(s') must have the same combinatorial structure as vr = svi{s) from s to Ui(s), 
it holds that vr' = s'ui(s') for any s' G a. EQ 

Lemmashows that the visibility for the vertex Vi{s) is preserved within cell a. Hence, we can 
simply say that Vi{s) is always visible from a or never visible from a. Since corners of V are also 
vertices of SPM(s), they are also always or never visible from a. 

Lemma 3. For any cell a of Aspm and z G {1, 2,... ,m}, if vertex Vi is visible from a, then Vi is 
a corner v G V ofV. 

Proof. If Vi is always visible from cr, then the shortest path from s to Vi{A is just the straight line 
segment svi{s) and therefore is unique. However, as discussed in Section^ the only way in which 
a vertex of SPM(s) has a single shortest path from s is when it is a corner of V. 0 


Lemma 4. For any cell a of Aspm and i G {1,2,..., m}, it holds that 

fi{s) = 

where Ui,Wi G V are two corners ofV uniquely determined by i. 


||s — Uj(s)|| if Vi is visible from a, 

len„._i„.(s,Uj(s)) otherwise, 
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1: Algorithm GeodesicCenter(P) 

2: Compute the SPM-equivalence decomposition Aspm of V. 

3: for each cell a of As pm do 

4: Specify the combinatorial structure of the shortest path maps SPM(s) for s € a. 

5: Identify the parameterized equations of the vertices of SPM(s). 

6: Let ui(s),..., Vm(s) be the parameterized points identified by the above step. 

7: Let fi(s) := d(s,Vi(s)) be the m bivariate functions for s G a. 

8: Compute the upper envelope of the m graphs {z = fi{s)}. 

9: Find all points with the lowest z-coordinate in Uu. 

10: Store them as the cr-constrained geodesic centers with its z-value. 

11: end for 

12: return All Co-’s having the smallest z-value as cen('P), and its z-value as rad('P). 

13: end Algorithm 

Figure 3: An 0(n^^"’“’^)-time algorithm for computing cen('P) and rad(P) of a polygonal domain V. 

Proof. The case in which Vi is visible follows from Lemma Thus, it suffices to consider the 
opposite case. Pick any point sq G cr, and consider a shortest path vr from sq to Ui(so). Let 
Ui G V and Wi G V he the first and the last corners of V along vr. Since Vi is not visible from a 
(and in particular from sq), no shortest path from sq to Uj(so) can be the straight line segment 
soVi{so). Therefore, such corners Ui and Wi must exist. This implies that /i(so) = d(s 0 ) 4'i(so)) = 
len^.^^.(so, Ui(so)). By the definition of the SPM-equivalence decomposition Aspm, for any s G a, 
the shortest paths from any s G a to Vi{s) have the same combinatorial structure. Therefore, the 
path whose first corner is Ui and last corner is Wi must also be a shortest path. Hence, the lemma 
follows. EQ 

By combining these two observations, we can explicitly construct the functions /i, / 2 ,..., /m, 
and exploit them to compute the geodesic radius and center. The pseudocode of our algorithm can 
be found in Figure 

Theorem 1. The algorithm described in Figure^ correctly computes the geodesic radius and center 
of a polygonal domain with n corners in time for any e > 0. 

Proof. The correctness follows from the discussion above. That is, any center of V corresponds to 
a minimum of the upper envelope of the functions fi, f 2 ,..., fm- 

In order to show the time bound, we need an efficient tool to compute the upper envelope of 
functions. Civen a collection of N algebraic surface patches in we can compute their lower (or 
upper) envelope in time using the algorithms of Halperin and Sharir (for d = 3) 

or of Sharir (for d > 3). Note that the complexity of the resulting envelope is bounded by 
0(Arrf-i+^). 

Recall that the coordinates of each vertex Vi{s) of SPM(s) is an algebraic function [^. Lemma 
implies the functions fi are algebraic, too. Thus, we can apply the above algorithms to compute 
the upper envelope of the graphs of fi, and obtain an explicit expression of <I>. 

In our case, we have N = m = 0{n), since any shortest path map SPM(s) has 0{n) complexity. 
Each function fi has two degrees of freedom (i.e., the coordinates of s within a), so the graph of fi 
lies in three-dimensional space. That is, the upper envelope Ida of the functions fi can be computed 
in ©(n^"’"'^) for any positive e. Once the envelope is computed, we can hnd the points with the 





lowest z-coordinate in IA„ in the same time bound by traversing all faces of the envelope 1/^. Any 
point that minimizes lAa is a candidate for a geodesic center, and its image will be its corresponding 
radius. 

Consequently, we spend time per cell a of .4 spm- Since ^spm consists of cells, 

we obtain the claimed time bound 0 

5 Concluding Remarks 

We have presented the first algorithm that computes the geodesic radius and center of a general 
polygonal domain with holes. The running time of our algorithm is large, but still comparable with 
those for other related problems. A bottleneck of our algorithm is to compute the SPM-equivalence 
decomposition ^spMj which is very complicated and not well understood. The best known upper 
bound on the complexity of .4 spm is 0(n^°), and it is known how to construct a polygonal domain 
whose decomposition ^spm has n(n^) complexity [^. Thus, a better analysis on the upper bound 
for ^SPM would directly lead to an improvement to our algorithm. 

Another approach for improvement would be to use a coarser subdivision, such as the SPT- 
equivalence decomposition [^. The SPT-equivalence decomposition .4 spt only requires shortest 
path trees SPT(s) for all s in each cell of ^spt to be isomorphic, rather than equivalence between 
SPM(s). The complexity of this subdivision is O(n^), which is much smaller than that of ^spM) 
and has similar (albeit slightly weaker) properties to those of »4 spm- Ideally, we would want an 
algorithm that can compute the cr-constrained geodesic radius for a cell a of ^spt in 
time so that overall the running time improves Theorem However, all of our attempts needed 
significantly more than Q{n^) time, which lead to even slower algorithms. 

Throughout the paper, we have focused on the exact computation of the geodesic radius and 
centers, but one could also consider the approximation variant. By the triangular inequality, 
any point s G "P and its farthest point t G V give a 2-approximation of the radius. That is, 
rad(P) < <h(s) < 2rad(P) for any s G V. Similarly, we can obtain a (1 -|- e)-approximation by 
using a standard grid technique: Scale V so that V fits into a unit square, and partition V with a 
grid of size 0{e~^) x 0(e“^). Define the set D to be the point set containing the grid points that 
are inside P, and intersection points between boundary edges and grid edges. We then observe 
that the distance between any two points s and t in P is within a (1 -|- e)-factor of the distance 
between two points of D that are closest from s and t, respectively. Hence, we conclude that 
min^go <I)(z) < (1 + e) •rad(P). Since D consists of 0{e~^{e~^ -\- n)) points and it takes 0(n log re) 
time per each point z G D to find its farthest neighbors using the shortest path map SPM( 2 ;), this 
algorithm runs in 0((^ + ^) log re) time0 No subquadratic-time approximation algorithm with 
factor less than 2 is known so far. 
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